
    ##############################################################################################
    # R code to replicate
    # Figure 2 of Measuring Foreign Policy Positions of Members of the US Congress
    ### Note: Run this file only after running Figure3_R_code.R first.
    ### You need "Summary_Outputs.csv", which will be created by MCMC estimation for Figure 3
    ##############################################################################################


    T <- 33  # number of congresses
    year.odd <- seq(from=1945, to=2010, by=2)

    Outputs <- read.csv("Summary_Outputs.csv", header=T) # Summary of Estimates from Figure 3 Estimation
    dim(Outputs)
    names(Outputs)
    
    # Storage    
    class_rates <- matrix(NA, nrow=T, ncol=4)

    ##############################################################
    ### Loop over congresses

    for (t in 1:T){

    # Merge Estimates with Roll Call Date    
    Vote_Data <- read.table(paste("All_FP_votes_s", t+78, "th.txt", sep=""), header=T)
    
    Measures <- subset(Outputs, time==t, select=c(ICPSR_ID, time, congress, IRT_mean, DW_NOM1, DW_NOM2))
        
    Data <- merge(Measures, Vote_Data, by="ICPSR_ID", all=F)
    dim(Data)
    
    Votes <- Data[ , 11:ncol(Data)]
    M <- ncol(Votes)
    N <- nrow(Votes)

    # Create GOP dummy and other IVs
    gop <- ifelse(Data$party==200, 1, 0)
    irt <- Data$IRT_mean
    d1n <- Data$DW_NOM1
    d2n <- Data$DW_NOM2


    #############################
    ### Compute classification rates for logit models with IRT measure as the only IV
    
    logit_class_irt <- rep(NA, M)
    
    for (j in 1:M){
    Y <- Votes[,j] 
    coeffs <- round(summary(glm(Y ~ irt, family='binomial'))$coeff, digits=2)   

    # (1) Classification Rate
    P <- rep(NA, N)
    for (i in 1:N){  # Loop over legislators
            P[i] <- plogis(coeffs[1,1] + coeffs[2,1]*irt[i])
    }

    ###### Generate Y.pred
    P[is.na(P)] <- 9 
    Y.pred <- rep(NA, N)
    for (i in 1:N){
    if (P[i] == 9) Y.pred[i] <- 9
    else {if (P[i] >= .5) Y.pred[i] <- 1
        else Y.pred[i] <- 0}
    }

    ### Classification Rates
    # Clean up data
    Y1 <- as.numeric(Y)
    Y1[is.na(Y1)] <- 9 # To change "NA" to "9". 

    # 1. For all senators
    Count <- rep(NA, N)
    for (i in 1:N){
    if (Y1[i] == 9) {Count[i] <- 9}
    else {if (Y1[i] == Y.pred[i]) {Count[i] <- 1}
          else {if (Y.pred[i] == 9) {Count[i] <- 9}
              else Count[i] <- 0}}
    }

    Count1 <- Count
    Count1[Count1==9] <- NA

    # to calculate the total number of 1's (=right classification)
    one.count <- rep(0, N)
    one.count[Count==1] <- 1

    # to calculate the total number of 0's (=wrong classification)
    zero.count <- rep(0, N)
    zero.count[Count==0] <- 1

    # Classification rate for all
    logit_class_irt[j] <- 100*(sum(one.count)/(sum(one.count)+ sum(zero.count)))
    } 
    
    class_rates[t,1] <- mean(logit_class_irt)
       

    #################################
    ### Compute classification rates for logit models with Dim 1 DW-NOM as the only IV
    
    logit_class_d1n <- rep(NA, M)
    
    for (j in 1:M){
    Y <- Votes[,j] 
    coeffs <- round(summary(glm(Y ~ d1n, family='binomial'))$coeff, digits=2)   

    # (1) Classification Rate
    P <- rep(NA, N)
    for (i in 1:N){  # Loop over legislators
            P[i] <- plogis(coeffs[1,1] + coeffs[2,1]*d1n[i])
    }
    ###### Generate Y.pred
    P[is.na(P)] <- 9 
    Y.pred <- rep(NA, N)
    for (i in 1:N){
    if (P[i] == 9) Y.pred[i] <- 9
    else {if (P[i] >= .5) Y.pred[i] <- 1
        else Y.pred[i] <- 0}
    }

    ### Classification Rates
    # Clean up data
    Y1 <- as.numeric(Y)
    Y1[is.na(Y1)] <- 9 # To change "NA" to "9". 
    
    # 1. For all senators
    Count <- rep(NA, N)
    for (i in 1:N){
    if (Y1[i] == 9) {Count[i] <- 9}
    else {if (Y1[i] == Y.pred[i]) {Count[i] <- 1}
          else {if (Y.pred[i] == 9) {Count[i] <- 9}
              else Count[i] <- 0}}
    }

    Count1 <- Count
    Count1[Count1==9] <- NA

    # to calculate the total number of 1's (=right classification)
    one.count <- rep(0, N)
    one.count[Count==1] <- 1
    sum(one.count)

    # to calculate the total number of 0's (=wrong classification)
    zero.count <- rep(0, N)
    zero.count[Count==0] <- 1
    sum(zero.count)

    # Classification rate for all
    logit_class_d1n[j] <- 100*(sum(one.count)/(sum(one.count)+ sum(zero.count)))
    } 

    class_rates[t,2] <- mean(logit_class_d1n)


    #################################
    ### Compute classification rates for logit models with Dim 2 DW-NOM as the only IV
    
    logit_class_d2n <- rep(NA, M)
    
    for (j in 1:M){
    Y <- Votes[,j] 
    coeffs <- round(summary(glm(Y ~ d2n, family='binomial'))$coeff, digits=2)   

    # (1) Classification Rate
    P <- rep(NA, N)
    for (i in 1:N){  # Loop over legislators
            P[i] <- plogis(coeffs[1,1] + coeffs[2,1]*d2n[i])
    }
    ###### Generate Y.pred
    P[is.na(P)] <- 9 
    Y.pred <- rep(NA, N)
    for (i in 1:N){
    if (P[i] == 9) Y.pred[i] <- 9
    else {if (P[i] >= .5) Y.pred[i] <- 1
        else Y.pred[i] <- 0}
    }

    ### Classification Rates
    # Clean up data
    Y1 <- as.numeric(Y)
    Y1[is.na(Y1)] <- 9 # To change "NA" to "9". 

    # 1. For all senators
    Count <- rep(NA, N)
    for (i in 1:N){
    if (Y1[i] == 9) {Count[i] <- 9}
    else {if (Y1[i] == Y.pred[i]) {Count[i] <- 1}
          else {if (Y.pred[i] == 9) {Count[i] <- 9}
              else Count[i] <- 0}}
    }

    Count1 <- Count
    Count1[Count1==9] <- NA

    # to calculate the total number of 1's (=right classification)
    one.count <- rep(0, N)
    one.count[Count==1] <- 1

    # to calculate the total number of 0's (=wrong classification)
    zero.count <- rep(0, N)
    zero.count[Count==0] <- 1

    # Classification rate for all
    logit_class_d2n[j] <- 100*(sum(one.count)/(sum(one.count)+ sum(zero.count)))
    } 

    class_rates[t,3] <- mean(logit_class_d2n)


    #################################
    ### Compute classification rates for logit models with Party dummy (GOP) as the only IV
    
    logit_class_gop <- rep(NA, M)
    
    for (j in 1:M){
    Y <- Votes[,j] 
    coeffs <- round(summary(glm(Y ~ gop, family='binomial'))$coeff, digits=2)   

    # (1) Classification Rate
    P <- rep(NA, N)
    for (i in 1:N){  # Loop over legislators
            P[i] <- plogis(coeffs[1,1] + coeffs[2,1]*gop[i])
    }
    ###### Generate Y.pred
    P[is.na(P)] <- 9 
    Y.pred <- rep(NA, N)
    for (i in 1:N){
    if (P[i] == 9) Y.pred[i] <- 9
    else {if (P[i] >= .5) Y.pred[i] <- 1
        else Y.pred[i] <- 0}
    }

    ### Classification Rates
    # Clean up data
    Y1 <- as.numeric(Y)
    Y1[is.na(Y1)] <- 9 # To change "NA" to "9". 

    # 1. For all senators
    Count <- rep(NA, N)
    for (i in 1:N){
    if (Y1[i] == 9) {Count[i] <- 9}
    else {if (Y1[i] == Y.pred[i]) {Count[i] <- 1}
          else {if (Y.pred[i] == 9) {Count[i] <- 9}
              else Count[i] <- 0}}
    }

    Count1 <- Count
    Count1[Count1==9] <- NA

    # to calculate the total number of 1's (=right classification)
    one.count <- rep(0, N)
    one.count[Count==1] <- 1
    sum(one.count)

    # to calculate the total number of 0's (=wrong classification)
    zero.count <- rep(0, N)
    zero.count[Count==0] <- 1
    sum(zero.count)

    # Classification rate for all
    logit_class_gop[j] <- 100*(sum(one.count)/(sum(one.count)+ sum(zero.count)))
    } 

    class_rates[t,4] <- mean(logit_class_gop)

    colnames(class_rates) <- c("IRT", "D1N", "D2N", "GOP")
    }
    #######################################################################
    ### the end of the loop over Congresses


    ### Present the results

    postscript("Figure2.eps", height=6, width=8, horizontal=F)
    par(mfrow=c(1, 1))
    plot(1:T, class_rates[,1], type='o', lty=1, col='blue', xaxt="n", xlab="", ylab="Overall Classification Success Rate (%)", ylim=c(50, 100))
    points(1:T, class_rates[,2],  type='o', lty=2, col='red')
    points(1:T, class_rates[,3],  type='o', lty=5, col='purple')
    points(1:T, class_rates[,4],  type='o', lty=4, col='black')
    axis(side=1, las=3, at=1:T, labels=year.odd, cex.axis=.8)
    legend('topleft', legend=c("Bayesian IRT", "DW-NOM Dim 1", "DW-NOM Dim 2", "Party Dummy"), lty=c(1, 2, 5, 4), col=c('blue', 'red', 'purple', 'black'),cex=.8)
    dev.off()
    

